Mock min and max temp list for testing
t_min_list = list("1"=-12, "2"=-8, "3"=-4, "4"=-3, "5"=8, "6"=10, "7"=11, "8"=15, "9"=18, "10"=3, "11"=-2, "12"=-10)
t_max_list = list("1"=-1, "2"=-3, "3"=-2, "4"=0, "5"=18, "6"=22, "7"=28, "8"=23, "9"=22, "10"=6, "11"=3, "12"=-7)
library(ggplot2)
library(viridis)
## Loading required package: viridisLite
library(ggExtra)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
NFFD
nffd_param <- read.csv(file = "../optimizedParameterTables/param_NFFD.csv", sep=',', header = TRUE)
# nffd: number of frost free days
# m: month of the year
# tm: min temperature for that month
nffd <- function(m, tm) {
optimized_params <- nffd_param[nffd_param$Month == m,]
a <- optimized_params$a
b <- optimized_params$b
t0 <- optimized_params$T0
cat(t0)
return( a/(1 + exp(-(tm - t0)/b)))
}
nffd(2, -3)
## 3.8
## [1] 0.03910693
Fucntion to compute NFFD with data
compute_nffd <- function (month, scenario, model) {
month_two_digits = sprintf("%02d",month)
file = paste("../data/ensmax.COM.Tmin", month_two_digits, ".csv", sep="")
tmin_month <- read.csv(file, sep=',', header = TRUE)
tmin_month <- tmin_month[tmin_month$scenario == 'historical',][c("Year",model)]
nffd_value <- lapply(tmin_month[model], function(tm) nffd(month, tm))
nffd_year <- tmin_month$Year
#nffd_date <- as.Date(with(tmin_month,paste(Year, month, "01", sep="-")),"%Y-%m-%d")
nffd_df <- data.frame(nffd_year,nffd_value)
names(nffd_df) <- c("Year", "Value")
return(nffd_df)
}
Test Monthly NFFD
scenario = "historical"
model = "ACCESS.ESM1.5"
jan_nffd <- compute_nffd(1, scenario, model)
## 3.82
feb_nffd <- compute_nffd(2, scenario, model)
## 3.8
mar_nffd <- compute_nffd(3, scenario, model)
## 3.6
apr_nffd <- compute_nffd(4, scenario, model)
## 3.13
may_nffd <- compute_nffd(5, scenario, model)
## 2.78
jun_nffd <- compute_nffd(6, scenario, model)
## 2.24
jul_nffd <- compute_nffd(7, scenario, model)
## 1.23
aug_nffd <- compute_nffd(8, scenario, model)
## 1.71
sep_nffd <- compute_nffd(9, scenario, model)
## 2.72
oct_nffd <- compute_nffd(10, scenario, model)
## 3.23
nov_nffd <- compute_nffd(11, scenario, model)
## 3.47
dec_nffd <- compute_nffd(12, scenario, model)
## 3.63
createTimeSeries(jan_nffd, "NFFD Historical Jan")
createTimeSeries(feb_nffd, "NFFD Historical Feb")
createTimeSeries(mar_nffd, "NFFD Historical Mar")
createTimeSeries(apr_nffd, "NFFD Historical Apr")
createTimeSeries(may_nffd, "NFFD Historical May")
createTimeSeries(jun_nffd, "NFFD Historical Jun")
createTimeSeries(jul_nffd, "NFFD Historical Jul")
createTimeSeries(aug_nffd, "NFFD Historical Aug")
createTimeSeries(sep_nffd, "NFFD Historical Sep")
createTimeSeries(oct_nffd, "NFFD Historical Oct")
createTimeSeries(nov_nffd, "NFFD Historical Nov")
createTimeSeries(dec_nffd, "NFFD Historical Dec")
Function to compute seasonal NFFD
#Winter
season_nffd <- function(months, scenario, model) {
list_of_monthly_nffd = list()
for(month in months) {
list_of_monthly_nffd[[month]] <- compute_nffd(month, scenario, model)
}
nffd_all <- data.table::rbindlist(list_of_monthly_nffd)
agg_summ_nffd_all = aggregate(Value~Year,nffd_all,FUN = sum)
return(agg_summ_nffd_all)
}
Compute seasonal NFFD’s
scenario = "historical"
model = "ACCESS.ESM1.5"
winter_months = c(12,1,2)
nffd_winter <- season_nffd(winter_months,scenario, model)
## 3.633.823.8
spring_months = c(3,4,5)
nffd_spring <- season_nffd(spring_months,scenario, model)
## 3.63.132.78
summer_months = c(6,7,8)
nffd_summer <- season_nffd(summer_months,scenario, model)
## 2.241.231.71
autumn_months = c(9,10,11)
nffd_autum <- season_nffd(autumn_months,scenario, model)
## 2.723.233.47
createTimeSeries(nffd_winter, "NFFD Historical Winter")
createTimeSeries(nffd_spring, "NFFD Historical Spring")
createTimeSeries(nffd_summer, "NFFD Historical Summer")
createTimeSeries(nffd_autum, "NFFD Historical Autum")
DD
dd_param_below_0 <- read.csv(file = "../optimizedParameterTables/param_DD_S1.csv", sep=',', header = TRUE)
dd_param_above_5 <- read.csv(file = "../optimizedParameterTables/param_DD_S2.csv", sep=',', header = TRUE)
dd_param_below_18 <- read.csv(file = "../optimizedParameterTables/param_DD_S3.csv", sep=',', header = TRUE)
dd_param_above_18 <- read.csv(file = "../optimizedParameterTables/param_DD_S4.csv", sep=',', header = TRUE)
dd <- function(m, tm) {
# ***********************************************
# Question, I think it should be dd < 5 (not > 5)
# ***********************************************
dd_param <- ''
if(tm < 0) {
dd_param <- dd_param_below_0
} else if(tm < 5) {
dd_param <- dd_param_above_5[dd_param_above_5$Region == "All"]
} else if(tm < 18) {
dd_param <- dd_param_below_18
} else {
dd_param <- dd_param_above_18
}
optimized_params <- dd_param[dd_param$Month == m,]
k <- optimized_params$a
a <- optimized_params$a
b <- optimized_params$b
t0 <- optimized_params$T0
c <- optimized_params$a
beta <- optimized_params$a
if(tm > k) {
return( a/(1 + exp(-(tm - t0)/b)))
} else {
return(c + (beta * tm))
}
}
dd(4, 14)
## [1] 4878.885
Test functions